genenral notes or findings:
- All series were found to be I(1), i.e. stationary after first differencing.
SetupΒΆ
# relevent libraries
# install.packages("quantmod")
# install.packages("fredr")
# install.packages("ggfortify")
# install.packages('urca')
# install.packages("tseries")
# install.packages("forecast")
# install.packages("dynlm")
# install.packages("stargazer")
# install.packages("pracma")
# install.packages("dLagM")
# install.packages("gets")
# install.packages("car")
# install.packages("lmtest")
# install.packages("vars")
# install.packages("tseries")
# install.packages("strucchange")
# install.packages("graphics")
# install.packages("grDevices")
# install.packages("tsDyn")
options(warn=-1)
library(xts)
library(zoo)
library(ggplot2)
library(ggfortify)
library(urca)
library(forecast)
data <- read.csv('dataset.csv')
head(data)
tail(data)
| Date | Europe.Brent.Spot.Price.FOB | Monthly.OECD.petroleum.and.other.liquids.stocks | World.Industrial.Production.index | OPEC.Production..million.b.d. | Baltic.Dry.Index..BADI. | |
|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | |
| 1 | Jan-00 | 25.51 | 3768.491 | NA | 24.18 | 1319 |
| 2 | Feb-00 | 27.78 | 3744.295 | NA | 24.95 | 1531 |
| 3 | Mar-00 | 27.49 | 3731.368 | NA | 24.91 | 1660 |
| 4 | Apr-00 | 22.76 | 3770.587 | NA | 25.69 | 1628 |
| 5 | May-00 | 27.74 | 3789.751 | NA | 26.25 | 1566 |
| 6 | Jun-00 | 29.80 | 3820.624 | NA | 25.83 | 1616 |
| Date | Europe.Brent.Spot.Price.FOB | Monthly.OECD.petroleum.and.other.liquids.stocks | World.Industrial.Production.index | OPEC.Production..million.b.d. | Baltic.Dry.Index..BADI. | |
|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | |
| 295 | Jul-24 | 85.15 | 4032.39 | 107.26 | 29.04 | 1708 |
| 296 | Aug-24 | 80.36 | 4047.26 | 107.38 | 28.94 | 1814 |
| 297 | Sep-24 | 74.02 | 4016.63 | 107.61 | 28.23 | 2084 |
| 298 | Oct-24 | 75.63 | 3981.54 | 108.16 | 28.71 | 1388 |
| 299 | Nov-24 | 74.35 | 3979.63 | 108.16 | 28.74 | 1354 |
| 300 | Dec-24 | 73.86 | 3975.76 | 109.04 | 28.82 | 997 |
# convert "Jan-90" to yearmon (monthly format)
d1 <- as.yearmon(data[,1], format = "%b-%y")
# convert to xts
Brent <- xts(data[,2], order.by = d1) # date-indexed series (Crude Oil Prices: Brent - Europe)
# we have no missing values on Brent, but we omit just in case
Brent = na.omit(Brent)
# we repeat the same process for the other series
Oil_Stocks <- xts(data[,3], order.by = d1)
Oil_Stocks = na.omit(Oil_Stocks)
WIP <- xts(data[,4], order.by = d1)
WIP = na.omit(WIP) # We do not have values of WIP before Jan 2016, all values before are omitted
OPEC <- xts(data[,5], order.by = d1)
OPEC = na.omit(OPEC)
BDI <- xts(data[,6], order.by = d1)
BDI = na.omit(BDI)
head(Brent)
tail(Brent)
[,1] Jan 2000 25.51 Feb 2000 27.78 Mar 2000 27.49 Apr 2000 22.76 May 2000 27.74 Jun 2000 29.80
[,1] Jul 2024 85.15 Aug 2024 80.36 Sep 2024 74.02 Oct 2024 75.63 Nov 2024 74.35 Dec 2024 73.86
Data goes from Jan 2000 to Dec 2024
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(Brent, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Europe Brent Spot Price FOB (Dollars per Barrel)")
fig
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(Oil_Stocks, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Monthly OECD petroleum and other liquids stocks - million barrels")
fig
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(WIP, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "World Industrial Production index")
fig
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(OPEC, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "OPEC Production (million b/d)")
fig
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(BDI, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Baltic Dry Index (BADI)")
fig
all_series <- merge(Brent, Oil_Stocks, WIP, OPEC, BDI)
colnames(all_series) <- c("Brent", "Oil Stocks", "World Industrial Production", "OPEC Production", "Baltic Dry Index")
df_long <- fortify(all_series, melt = TRUE)
colnames(df_long) <- c("Date", "Series", "Value")
options(repr.plot.width = 30, repr.plot.height = 20)
fig <- ggplot(df_long, aes(x = Date, y = Value, colour = Series)) +
geom_line(size = 1) +
theme_light() +
theme(
# aspect.ratio = 1,
plot.margin = ggplot2::margin(0.2, 0.2, 0.2, 0.2, "cm"),
text = element_text(size = 30),
legend.position = c(0.15, 0.85),
legend.background = element_rect(fill = alpha("white", 0.6)),
legend.key = element_rect(fill = NA)
) +
labs(
x = "Month - Year",
y = "Value",
colour = "Series"
)
fig
ADF Testing of Brent Oil PricesΒΆ
summary(ur.df(Brent, type='trend', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-21.8417 -3.9894 0.6469 3.7461 17.4151
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.902228 1.170712 2.479 0.01380 *
z.lag.1 -0.041953 0.014119 -2.971 0.00324 **
tt 0.001226 0.004833 0.254 0.80003
z.diff.lag 0.334422 0.057681 5.798 1.91e-08 ***
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 5.948 on 265 degrees of freedom
Multiple R-squared: 0.1304, Adjusted R-squared: 0.1206
F-statistic: 13.25 on 3 and 265 DF, p-value: 4.395e-08
Value of test-statistic is: -2.9714 3.0686 4.5485
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.98 -3.42 -3.13
phi2 6.15 4.71 4.05
phi3 8.34 6.30 5.36
We cannot reject the null hypothesis that the constant and deterministic trend are jointly zero because $\varphi_2=3.07<4.71$.
optional: (We also cannot reject that the deterministic trend is zero because $\varphi_3=4.55<6.30$.)
Therefore we re-run the ADF test with specification = "none" (no constant, no determistic trend).
summary(ur.df(Brent, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-20.894 -3.198 1.357 4.168 17.064
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.003547 0.004836 -0.733 0.464
z.diff.lag 0.320206 0.058071 5.514 8.27e-08 ***
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 6.022 on 267 degrees of freedom
Multiple R-squared: 0.1027, Adjusted R-squared: 0.09599
F-statistic: 15.28 on 2 and 267 DF, p-value: 5.204e-07
Value of test-statistic is: -0.7334
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We cannot reject the null that there exists a unit root at the 5pct significance level since $\tau_1=-0.7334>-1.095"$.
Therefore we difference the series. Chatgpt why we also take the log:
We log the Brent crude oil price series because crude oil prices exhibit strong multiplicative volatility: they move in percentages rather than fixed amounts, and large level increases create proportionally larger fluctuations. Taking the logarithm stabilizes the variance, making the series more homoskedastic and better aligned with the linear assumptions of ARIMA and ADL models. Log-transforming also converts differences into approximate percentage changes, so the series becomes interpretable as monthly returns, which are typically stationary even when price levels are not. This improves model performance, reduces the influence of extreme outliers, and produces more reliable forecasts compared with using raw price levels.
- Oil prices are strictly positive and highly volatile β logs stabilize variance
Brent has:
multiplicative shocks (prices jump by percentages, not by fixed amounts)
strong heteroskedasticity
explosive episodes (2008, COVID crash, 2022 spike)
Taking logs makes the series behave more like a linear, additive process, which is exactly what ARIMA/ADL/VAR models assume.
Without logging:
large spikes dominate the model
residuals become heteroskedastic
parameters can be unstable
forecasts become scale-dependent
This is all bad.
Brent_log <- log(Brent)
dBrent_log = diff(Brent_log)
# diff() always produces an NA in the first observation (there is no lagged value for t=1).
dBrent_log <- na.omit(dBrent_log)
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(dBrent_log, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Transformed Europe Brent Spot Price FOB (Dollars per Barrel)")
fig
summary(ur.df(dBrent_log, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.49795 -0.05073 0.01671 0.06067 0.54964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.96292 0.13098 -7.352 2.5e-12 ***
z.diff.lag1 0.26852 0.11993 2.239 0.0260 *
z.diff.lag2 0.13147 0.10510 1.251 0.2121
z.diff.lag3 0.08946 0.09048 0.989 0.3237
z.diff.lag4 -0.01551 0.07460 -0.208 0.8355
z.diff.lag5 0.11071 0.06115 1.810 0.0714 .
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.0995 on 262 degrees of freedom
Multiple R-squared: 0.3991, Adjusted R-squared: 0.3854
F-statistic: 29.01 on 6 and 262 DF, p-value: < 2.2e-16
Value of test-statistic is: -7.3519
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We reject the null of a unit root, meaning that the series is stationary. I.e. there is enough evidence to conclude that the series is stationary because $\tau_1=-7.35<-1.95$.
ARIMA Subset Search with AIC/BICΒΆ
p = 4 # Maximum AR order
q = 4 # Maximum MA order
pq = (p+1) * (q+1) # Total number of (p,q) combinations including 0
resAIC = rep(0, pq) # Storage for AIC values
resBIC = rep(0, pq) # Storage for BIC values
pdf = rep(0, pq) # Storage for AR order (p)
qdf = rep(0, pq) # Storage for MA order (q)
# if we condisider subarima models, where coefficients can be or not be, then we jump to 2^(p+q) models. but this function does not consider this. the next section does.
k = 0 # Counter for results
for (i in 0:p) {
for (j in 0:q) {
k = k + 1
# Estimate ARMA(p,q) model for stationary series dy
model = arima(dBrent_log, order = c(i,0,j), include.mean = FALSE)
# Store information criteria
resAIC[k] = AIC(model)
resBIC[k] = BIC(model)
# Store corresponding lag orders
pdf[k] = i
qdf[k] = j
}
}
# Combine and display results: p, q, AIC, and BIC
results = cbind(pdf, qdf, data.frame(resAIC), data.frame(resBIC))
results
# Find the model with the lowest AIC
bestAIC = results[which.min(results$resAIC), ] # Row with min AIC
# Find the model with the lowest BIC
bestBIC = results[which.min(results$resBIC), ] # Row with min BIC
# Display them
print("According to the AIC information criterion, the best model is")
bestAIC
print("According to the BIC information criterion, the best model is")
bestBIC
| qdf | resAIC | resBIC | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> |
| 0 | 0 | -501.7018 | -498.0013 |
| 0 | 1 | -520.8478 | -513.4469 |
| 0 | 2 | -520.4051 | -509.3038 |
| 0 | 3 | -518.5640 | -503.7622 |
| 0 | 4 | -521.7599 | -503.2576 |
| 1 | 0 | -516.5980 | -509.1971 |
| 1 | 1 | -520.5651 | -509.4637 |
| 1 | 2 | -521.0218 | -506.2201 |
| 1 | 3 | -519.1545 | -500.6523 |
| 1 | 4 | -520.5130 | -498.3103 |
| 2 | 0 | -520.7810 | -509.6797 |
| 2 | 1 | -519.8568 | -505.0550 |
| 2 | 2 | -519.2473 | -500.7450 |
| 2 | 3 | -520.8788 | -498.6762 |
| 2 | 4 | -518.5790 | -492.6759 |
| 3 | 0 | -519.1719 | -504.3701 |
| 3 | 1 | -520.7860 | -502.2838 |
| 3 | 2 | -520.6640 | -498.4614 |
| 3 | 3 | -516.4526 | -490.5495 |
| 3 | 4 | -531.3006 | -501.6970 |
| 4 | 0 | -519.1070 | -500.6047 |
| 4 | 1 | -520.1893 | -497.9866 |
| 4 | 2 | -531.5167 | -505.6136 |
| 4 | 3 | -528.8733 | -499.2698 |
| 4 | 4 | -527.6806 | -494.3767 |
[1] "According to the AIC information criterion, the best model is"
| qdf | resAIC | resBIC | ||
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | |
| 23 | 4 | 2 | -531.5167 | -505.6136 |
[1] "According to the BIC information criterion, the best model is"
| qdf | resAIC | resBIC | ||
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | |
| 2 | 0 | 1 | -520.8478 | -513.4469 |
fit_uni <- auto.arima(dBrent_log,
max.d = 0, # no differencing allowed (series already made stationary)
max.p = 4, # maximum AR order (p)
max.q = 4, # maximum MA order (q)
max.D = 0, # no seasonal differencing
max.P = 0, # no seasonal AR terms
max.Q = 0, # no seasonal MA terms
ic = "aic", # choose the model that minimizes AIC
stepwise = FALSE, # search over all models (not just stepwise path)
approximation = FALSE # use exact maximum likelihood (slower but more accurate)
)
summary(fit_uni)
Series: dBrent_log
ARIMA(0,0,4) with zero mean
Coefficients:
ma1 ma2 ma3 ma4
0.2754 -0.0875 -0.0653 -0.1317
s.e. 0.0579 0.0591 0.0631 0.0576
sigma^2 = 0.01002: log likelihood = 265.88
AIC=-521.76 AICc=-521.56 BIC=-503.26
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.003633379 0.09941295 0.07207996 94.73877 169.993 0.692788
ACF1
Training set -0.005245346
We evaluated autoregressive moving average (ARMA) models with up to 4 lags in both the autoregressive and moving average terms. This meant that a total of 25 models (all combinations of p = 0,β¦,4 and q = 0,β¦,4) were estimated. The choice of 4 lags reflects a balance between allowing enough dynamics to capture quarterly fluctuations in GDP and keeping the model space manageable for interpretation.
The models were compared using the Akaike Information Criterion (AIC), and the search considered all possible specifications within this range (an exact search rather than a stepwise approximation). Likelihood estimation was carried out without approximation. Based on this procedure, the ARMA(0,0) model, which corresponds to white noise, was selected as the best fitting model.
INCLUDE THIS:
An exhaustive ARMA(p,q) search with AIC favored a higher-order ARMA(4,2) model, which is common when using AIC because its penalty for additional parameters is relatively mild; however, auto.arima, which uses AICc and enforces stationarity/invertibility constraints, selected a more parsimonious ARMA(0,4) model. We therefore use the auto.arima model as the benchmark, consistent with forecasting practice and the course methodology.
Finally, we decided to let auto.arima find the best model without restrictions so we ran:
fit_uni <- auto.arima(dBrent_log,
max.d = 0, # no differencing allowed (series already made stationary)
# max.p = 4, # maximum AR order (p)
# max.q = 4, # maximum MA order (q)
# max.D = 0, # no seasonal differencing
# max.P = 0, # no seasonal AR terms
# max.Q = 0, # no seasonal MA terms
ic = "aic", # choose the model that minimizes AIC
stepwise = FALSE, # search over all models (not just stepwise path)
approximation = FALSE # use exact maximum likelihood (slower but more accurate)
)
summary(fit_uni)
Series: dBrent_log
ARIMA(0,0,4) with zero mean
Coefficients:
ma1 ma2 ma3 ma4
0.2754 -0.0875 -0.0653 -0.1317
s.e. 0.0579 0.0591 0.0631 0.0576
sigma^2 = 0.01002: log likelihood = 265.88
AIC=-521.76 AICc=-521.56 BIC=-503.26
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.003633379 0.09941295 0.07207996 94.73877 169.993 0.692788
ACF1
Training set -0.005245346
Interpretation of the Univariate Benchmark ModelΒΆ
Using exact maximum likelihood (approximation = FALSE), auto.arima selects an
$$ \text{ARIMA}(0,0,4) $$
model for the stationary series of log-differenced Brent prices,
$$ dBrent\_{log,t} = \Delta \log(Brent_t), $$
with zero mean. Because the series has already been differenced and is centered
around zero, the model contains only MA terms, capturing short-run, shock-driven
dynamics typical of commodity returns.
The estimated model is:
$$ dBrent\_{log,t} = \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \theta_3 \varepsilon_{t-3} + \theta_4 \varepsilon_{t-4} + \varepsilon_t, $$
with the following coefficients:
- $\theta_1 = 0.2754$ (significant, positive short-run persistence)
- $\theta_2 = -0.0875$
- $\theta_3 = -0.0653$
- $\theta_4 = -0.1317$
These MA terms reflect how past shocks (news surprises) propagate through monthly oil returns. The absence of AR terms indicates that the series does not exhibit meaningful autoregressive structure at monthly frequency once log-differenced.
The fit is strong:
- Residual variance: $$\sigma^2 \approx 0.0100$$
- AIC and AICc are very low (β521.76 and β521.56), indicating good in-sample fit
- BIC = β503.26 also supports the modelβs parsimony
- The first residual autocorrelation is essentially zero (ACF1 β β0.005), showing that the MA(4) structure captures all meaningful time dependence.
Overall, this ARIMA(0,0,4) model is a clean and parsimonious univariate benchmark: it treats monthly Brent log-returns as a short-memory process driven by the last few shocks, with no trend, drift, or autoregressive components. This serves as the baseline against which multivariate models (ADL or VAR with macroeconomic, financial, and supply-side predictors) will later be compared.
This ARIMA model forms our univariate benchmark, against which we later compare multivariate models incorporating macroeconomic, financial, and supply-side predictors.
# dBrent_log is your xts series of log-returns
dates <- index(dBrent_log)
# Fitted values from ARIMA, coerced to numeric and put on same dates
fitted_xts <- xts(as.numeric(fitted(fit_uni)), order.by = dates)
# Combine actual and fitted into one xts object
plot_xts <- merge(Actual = dBrent_log, Fitted = fitted_xts)
# Build data frame for ggplot
plot_df <- data.frame(
Date = index(plot_xts),
Actual = as.numeric(plot_xts$Actual),
Fitted = as.numeric(plot_xts$Fitted)
)
options(repr.plot.width = 20, repr.plot.height = 10)
fig <- ggplot(plot_df, aes(x = Date)) +
geom_line(aes(y = Actual, colour = "Actual"), linewidth = 0.7) +
geom_line(aes(y = Fitted, colour = "Fitted"), linewidth = 0.6) +
scale_colour_manual(values = c("Actual" = "black", "Fitted" = "red")) +
theme_light() +
theme(
plot.margin = ggplot2::margin(1, 1, 1, 1, "cm"),
text = element_text(size = 25),
legend.title = element_blank()
) +
labs(
title = "ARIMA(0,0,4) β Actual vs Fitted dBrent_log",
x = "Month - Year",
y = "dBrent_log"
)
fig
dim(dBrent_log)
dim(fitted_xts)
head(dBrent_log)
head(fitted_xts)
- 299
- 1
- 299
- 1
[,1] Feb 2000 0.08524581 Mar 2000 -0.01049404 Apr 2000 -0.18881769 May 2000 0.19787081 Jun 2000 0.07163298 Jul 2000 -0.03830838
[,1] Feb 2000 0.004154238 Mar 2000 0.019868432 Apr 2000 -0.017301407 May 2000 -0.050741019 Jun 2000 0.073614611 Jul 2000 -0.006555693
In an MA(4) model, the fitted value at time t depends on the previous four shocks, which are not directly observable at the start of the sample. However, auto.arima estimates these initial shocks using maximum likelihood and backcasting, allowing the model to produce fitted values for the entire sample without introducing missing observations. This is standard in ARIMA implementations: initial states are treated as parameters and inferred jointly with the MA coefficients. As a result, fitted(fit_uni) contains valid values from the first period onward, even though the early fitted points rely more heavily on the estimated pre-sample shocks.
Although the ARIMA model on log-returns cannot reproduce historical price levels, its forecasts remain fully meaningful: predicted log-returns can be accumulated and exponentiated to yield price forecasts. This is standard practice in financial econometrics, where stationary return models are used to generate price forecasts through the transformation (P_{t+h} = P_t \exp(\sum_{i=1}^h \widehat{d\log P}_{t+i})).
tail(Brent)
# price forecasts
fc <- forecast(fit_uni, h = 4)
last_price <- as.numeric(last(Brent))
fc_prices <- last_price * exp(cumsum(fc$mean))
fc_prices
[,1] Jul 2024 85.15 Aug 2024 80.36 Sep 2024 74.02 Oct 2024 75.63 Nov 2024 74.35 Dec 2024 73.86
- 74.5082607531298
- 74.2970033763458
- 74.6859194252936
- 74.7389284720804
ADF Testing of RegressorsΒΆ
We test the stationarity of the regressors and apply required transformations to them.
# Oil_Stocks, WIP, OPEC, BDI
Monthly OECD petroleum and other liquids stocksΒΆ
summary(ur.df(Oil_Stocks, type='trend', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-83.42 -15.80 0.74 17.53 133.90
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 98.438477 43.394678 2.268 0.024145 *
z.lag.1 -0.023730 0.010793 -2.199 0.028808 *
tt 0.011855 0.028938 0.410 0.682396
z.diff.lag1 0.157975 0.061095 2.586 0.010277 *
z.diff.lag2 0.149594 0.057928 2.582 0.010374 *
z.diff.lag3 0.202172 0.058366 3.464 0.000625 ***
z.diff.lag4 -0.037399 0.059724 -0.626 0.531747
z.diff.lag5 -0.076052 0.059790 -1.272 0.204543
z.diff.lag6 0.008252 0.059923 0.138 0.890579
z.diff.lag7 -0.103241 0.059110 -1.747 0.081921 .
z.diff.lag8 0.067807 0.059440 1.141 0.255045
z.diff.lag9 0.034403 0.059762 0.576 0.565345
z.diff.lag10 0.064062 0.059674 1.074 0.284055
z.diff.lag11 -0.037634 0.058673 -0.641 0.521830
z.diff.lag12 0.381223 0.058390 6.529 3.6e-10 ***
z.diff.lag13 -0.156074 0.062761 -2.487 0.013535 *
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 30.26 on 253 degrees of freedom
Multiple R-squared: 0.2983, Adjusted R-squared: 0.2567
F-statistic: 7.169 on 15 and 253 DF, p-value: 4.498e-13
Value of test-statistic is: -2.1987 1.8629 2.7943
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.98 -3.42 -3.13
phi2 6.15 4.71 4.05
phi3 8.34 6.30 5.36
We cannot reject the null hypothesis that the constant and deterministic trend are jointly zero because $\varphi_2=1.86<4.71$.
optional: (We also cannot reject that the deterministic trend is zero because $\varphi_3=2.79<6.30$.)
Therefore we re-run the ADF test with specification = "none" (no constant, no determistic trend).
summary(ur.df(Oil_Stocks, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-77.57 -17.42 1.55 16.02 128.44
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.0000493 0.0004387 -0.112 0.91061
z.diff.lag1 0.1592557 0.0614759 2.591 0.01013 *
z.diff.lag2 0.1391212 0.0578646 2.404 0.01692 *
z.diff.lag3 0.1915701 0.0583661 3.282 0.00117 **
z.diff.lag4 -0.0527497 0.0595435 -0.886 0.37650
z.diff.lag5 -0.0920165 0.0595558 -1.545 0.12358
z.diff.lag6 -0.0069859 0.0597194 -0.117 0.90697
z.diff.lag7 -0.1144304 0.0591365 -1.935 0.05409 .
z.diff.lag8 0.0586134 0.0595379 0.984 0.32582
z.diff.lag9 0.0238189 0.0598039 0.398 0.69075
z.diff.lag10 0.0541752 0.0597580 0.907 0.36549
z.diff.lag11 -0.0516730 0.0584254 -0.884 0.37730
z.diff.lag12 0.3651596 0.0579236 6.304 1.27e-09 ***
z.diff.lag13 -0.1833568 0.0615107 -2.981 0.00315 **
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 30.47 on 255 degrees of freedom
Multiple R-squared: 0.2828, Adjusted R-squared: 0.2434
F-statistic: 7.182 on 14 and 255 DF, p-value: 1.592e-12
Value of test-statistic is: -0.1124
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We cannot reject the null that there exists a unit root at the 5pct significance level since $\tau_1=-0.11>-1.95$.
Therefore we difference the series and take the log as the series has a clear upward trend and its variance grows with the level. Log-differencing stabilizes the variance and turns level changes into approximate percentage changes, making the transformed series much closer to weakly stationary. If you only difference levels, you often remove the trend but leave heteroskedasticity and scale effects, which leads to poorer ARIMA/VAR performance.
Oil_Stocks_log <- log(Oil_Stocks)
dOil_Stocks_log = diff(Oil_Stocks_log)
# diff() always produces an NA in the first observation (there is no lagged value for t=1).
dOil_Stocks_log <- na.omit(dOil_Stocks_log)
summary(ur.df(dOil_Stocks_log, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.0196585 -0.0039973 0.0002265 0.0037479 0.0282856
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.518077 0.151567 -3.418 0.000734 ***
z.diff.lag1 -0.344925 0.152969 -2.255 0.024990 *
z.diff.lag2 -0.203051 0.147825 -1.374 0.170777
z.diff.lag3 -0.008394 0.142398 -0.059 0.953042
z.diff.lag4 -0.063899 0.135330 -0.472 0.637207
z.diff.lag5 -0.157040 0.127390 -1.233 0.218806
z.diff.lag6 -0.166712 0.117649 -1.417 0.157696
z.diff.lag7 -0.280178 0.110117 -2.544 0.011538 *
z.diff.lag8 -0.228228 0.103512 -2.205 0.028358 *
z.diff.lag9 -0.205143 0.098066 -2.092 0.037439 *
z.diff.lag10 -0.144993 0.092872 -1.561 0.119712
z.diff.lag11 -0.189079 0.081131 -2.331 0.020557 *
z.diff.lag12 0.172815 0.061580 2.806 0.005397 **
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.007212 on 255 degrees of freedom
Multiple R-squared: 0.6023, Adjusted R-squared: 0.5821
F-statistic: 29.71 on 13 and 255 DF, p-value: < 2.2e-16
Value of test-statistic is: -3.4181
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
The transformed series is now stationary since we reject the null hypothesis of the presence of a unit root because $\tau_1=-3.42<-2.58$.
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(dOil_Stocks_log, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Transformed Monthly OECD petroleum and other liquids stocks")
fig
World Industrial Production index (WIP)ΒΆ
summary(ur.df(WIP, type='trend', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-8.9913 -0.2618 0.0944 0.5366 4.9502
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.85439 4.75024 3.338 0.001332 **
z.lag.1 -0.18306 0.05406 -3.386 0.001144 **
tt 0.03832 0.01305 2.938 0.004425 **
z.diff.lag 0.39480 0.10670 3.700 0.000415 ***
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 1.549 on 73 degrees of freedom
Multiple R-squared: 0.2204, Adjusted R-squared: 0.1884
F-statistic: 6.881 on 3 and 73 DF, p-value: 0.0003813
Value of test-statistic is: -3.3862 4.0205 5.8195
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
We cannot reject the null hypothesis that the constant and deterministic trend are jointly zero because $\varphi_2=4.02<4.75$.
optional: (We also cannot reject that the deterministic trend is zero because $\varphi_3=5.82<6.49$.)
Therefore we re-run the ADF test with specification = "none" (no constant, no determistic trend).
summary(ur.df(WIP, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-8.0803 -0.3883 0.0152 0.5767 6.8035
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 0.001254 0.001849 0.678 0.49993
z.diff.lag1 0.376357 0.113882 3.305 0.00147 **
z.diff.lag2 -0.212485 0.113944 -1.865 0.06617 .
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 1.62 on 74 degrees of freedom
Multiple R-squared: 0.1433, Adjusted R-squared: 0.1085
F-statistic: 4.125 on 3 and 74 DF, p-value: 0.009227
Value of test-statistic is: 0.6779
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We cannot reject the null that there exists a unit root at the 5pct significance level since $\tau_1=-0.68>-1.95$. The World Industrial Production Index displays a clear upward trend before and after the COVID-19 shock, with a large temporary collapse in early 2020. Because the level of the series grows over time and the variance appears to increase with the level, the process is unlikely to be stationary in levels. To obtain a stationary transformation suitable for ARIMA/ADL modelling, the series should be log-transformed and then differenced. Log-differencing stabilizes the variance and converts growth in the index into approximate percentage changes, producing a series that fluctuates around a stable mean without a deterministic trend.
WIP_log <- log(WIP)
dWIP_log = diff(WIP_log)
# diff() always produces an NA in the first observation (there is no lagged value for t=1).
dWIP_log <- na.omit(dWIP_log)
summary(ur.df(dWIP_log, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.0023731 -0.0004836 0.0000248 0.0004469 0.0033737
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.519391 0.152078 -3.415 0.000741 ***
z.diff.lag1 -0.346175 0.153434 -2.256 0.024907 *
z.diff.lag2 -0.204046 0.148273 -1.376 0.169982
z.diff.lag3 -0.009196 0.142813 -0.064 0.948709
z.diff.lag4 -0.065071 0.135697 -0.480 0.631971
z.diff.lag5 -0.158288 0.127708 -1.239 0.216317
z.diff.lag6 -0.168202 0.117925 -1.426 0.154993
z.diff.lag7 -0.281667 0.110352 -2.552 0.011281 *
z.diff.lag8 -0.230547 0.103720 -2.223 0.027107 *
z.diff.lag9 -0.207517 0.098251 -2.112 0.035649 *
z.diff.lag10 -0.146693 0.093030 -1.577 0.116074
z.diff.lag11 -0.189948 0.081234 -2.338 0.020146 *
z.diff.lag12 0.171498 0.061588 2.785 0.005761 **
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.0008643 on 255 degrees of freedom
Multiple R-squared: 0.6029, Adjusted R-squared: 0.5826
F-statistic: 29.78 on 13 and 255 DF, p-value: < 2.2e-16
Value of test-statistic is: -3.4153
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
The transformed series is now stationary since we reject the null hypothesis of the presence of a unit root because $\tau_1=-3.42<-1.95$.
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(dWIP_log, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Transformed World Industrial Production index")
fig
OPEC Production (million b/d)ΒΆ
summary(ur.df(OPEC, type='trend', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-5.7123 -0.2018 0.0145 0.2857 1.9667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.9282081 0.5657624 3.408 0.000756 ***
z.lag.1 -0.0620135 0.0192209 -3.226 0.001412 **
tt -0.0004843 0.0004496 -1.077 0.282385
z.diff.lag1 -0.0200903 0.0597006 -0.337 0.736748
z.diff.lag2 -0.1208657 0.0596477 -2.026 0.043737 *
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.5667 on 264 degrees of freedom
Multiple R-squared: 0.06198, Adjusted R-squared: 0.04777
F-statistic: 4.361 on 4 and 264 DF, p-value: 0.001972
Value of test-statistic is: -3.2264 4.212 6.1307
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.98 -3.42 -3.13
phi2 6.15 4.71 4.05
phi3 8.34 6.30 5.36
We cannot reject the null hypothesis that the constant and deterministic trend are jointly zero because $\varphi_2=4.21<4.71$.
optional: (We also cannot reject that the deterministic trend is zero because $\varphi_3=6.13<6.30$.)
Therefore we re-run the ADF test with specification = "none" (no constant, no determistic trend).
summary(ur.df(OPEC, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-5.7916 -0.1894 0.0262 0.2684 2.3263
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 0.0004771 0.0011952 0.399 0.6901
z.diff.lag1 -0.0315418 0.0603874 -0.522 0.6019
z.diff.lag2 -0.1315419 0.0603715 -2.179 0.0302 *
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.5778 on 266 degrees of freedom
Multiple R-squared: 0.01859, Adjusted R-squared: 0.00752
F-statistic: 1.679 on 3 and 266 DF, p-value: 0.1718
Value of test-statistic is: 0.3992
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We cannot reject the null that there exists a unit root at the 5pct significance level since $\tau_1=0.3992>-1.95$. OPEC crude oil production is clearly non-stationary in levels. The series exhibits long swings, persistent upward and downward movements, and clear structural breaks β most notably the collapse during the 2008β09 recession and the extreme COVID-19 drop in 2020. These level shifts imply that both the mean and the variance change over time, violating weak stationarity. To obtain a transformation suitable for ARIMA or ADL modelling, the series should be log-transformed and then differenced. Log-differencing stabilizes the scale, converts changes in production into approximate percentage changes, and produces a growth-rate series that fluctuates around a constant mean and is therefore much closer to stationarity.
OPEC_log <- log(OPEC)
dOPEC_log = diff(OPEC_log)
# diff() always produces an NA in the first observation (there is no lagged value for t=1).
dOPEC_log <- na.omit(dOPEC_log)
summary(ur.df(dOPEC_log, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.211770 -0.005906 0.001522 0.009675 0.093641
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -1.14594 0.08534 -13.428 <2e-16 ***
z.diff.lag 0.14847 0.05997 2.476 0.0139 *
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.02068 on 266 degrees of freedom
Multiple R-squared: 0.5101, Adjusted R-squared: 0.5065
F-statistic: 138.5 on 2 and 266 DF, p-value: < 2.2e-16
Value of test-statistic is: -13.4285
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
The transformed series is now stationary since we reject the null hypothesis of the presence of a unit root because $\tau_1=-13.43<-1.95$.
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(dOPEC_log, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Transformed OPEC Production (million b/d)")
fig
Baltic Dry Index (BADI)ΒΆ
summary(ur.df(BDI, type='trend', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression trend
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-2322.23 -235.44 -35.84 215.90 2770.72
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 320.662483 125.055729 2.564 0.010901 *
z.lag.1 -0.068877 0.022234 -3.098 0.002162 **
tt -0.976491 0.527081 -1.853 0.065059 .
z.diff.lag1 0.222114 0.058754 3.780 0.000194 ***
z.diff.lag2 0.001097 0.060312 0.018 0.985501
z.diff.lag3 0.088752 0.060353 1.471 0.142610
z.diff.lag4 -0.246625 0.059765 -4.127 4.95e-05 ***
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 595.3 on 262 degrees of freedom
Multiple R-squared: 0.1477, Adjusted R-squared: 0.1282
F-statistic: 7.567 on 6 and 262 DF, p-value: 1.691e-07
Value of test-statistic is: -3.0978 3.2749 4.9123
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.98 -3.42 -3.13
phi2 6.15 4.71 4.05
phi3 8.34 6.30 5.36
We cannot reject the null hypothesis that the constant and deterministic trend are jointly zero because $\varphi_2=3.27<4.71$.
optional: (We also cannot reject that the deterministic trend is zero because $\varphi_3=4.91<6.30$.)
Therefore we re-run the ADF test with specification = "none" (no constant, no determistic trend).
summary(ur.df(BDI, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-2479.85 -171.08 51.38 295.51 2641.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.01924 0.01242 -1.550 0.122359
z.diff.lag1 0.20508 0.05897 3.478 0.000592 ***
z.diff.lag2 -0.02245 0.06024 -0.373 0.709661
z.diff.lag3 0.06657 0.06038 1.103 0.271212
z.diff.lag4 -0.27434 0.05943 -4.616 6.11e-06 ***
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 601.3 on 264 degrees of freedom
Multiple R-squared: 0.1237, Adjusted R-squared: 0.1071
F-statistic: 7.454 on 5 and 264 DF, p-value: 1.466e-06
Value of test-statistic is: -1.5499
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
We cannot reject the null that there exists a unit root at the 5pct significance level since $\tau_1=-1.5499>-1.95$. The Baltic Dry Index (BDI) is strongly non-stationary in levels. The series displays large boomβbust cycles (notably the 2003β2008 surge and the collapse that followed), long spells of persistent drift, and clear breaks around major global demand shocks such as 2009 and 2020. Both the mean and the variance change substantially over time, violating the assumptions of weak stationarity. To model the growth rate of shipping activity in a meaningful way, the series should be log-transformed and then differenced. The log-first-difference transformation stabilizes the scale, turns level swings into approximate percentage changes, and produces a series that fluctuates around a stable mean, making it far more suitable for ARIMA, ADL, or VAR analysis.
BDI_log <- log(BDI)
dBDI_log = diff(BDI_log)
# diff() always produces an NA in the first observation (there is no lagged value for t=1).
dBDI_log <- na.omit(dBDI_log)
summary(ur.df(dBDI_log, type='none', lags=30, selectlags="AIC"))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-1.32804 -0.13536 0.02243 0.15470 1.18551
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -1.65427 0.20246 -8.171 1.35e-14 ***
z.diff.lag1 0.63720 0.18389 3.465 0.000620 ***
z.diff.lag2 0.56429 0.16546 3.410 0.000752 ***
z.diff.lag3 0.56560 0.14611 3.871 0.000137 ***
z.diff.lag4 0.36894 0.12696 2.906 0.003976 **
z.diff.lag5 0.29367 0.10946 2.683 0.007769 **
z.diff.lag6 0.22010 0.08751 2.515 0.012505 *
z.diff.lag7 0.14624 0.06183 2.365 0.018767 *
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Residual standard error: 0.2683 on 260 degrees of freedom
Multiple R-squared: 0.5276, Adjusted R-squared: 0.5131
F-statistic: 36.3 on 8 and 260 DF, p-value: < 2.2e-16
Value of test-statistic is: -8.1709
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
The transformed series is now stationary since we reject the null hypothesis of the presence of a unit root because $\tau_1=-8.1709<-1.95$.
options(repr.plot.width=30, repr.plot.height=20)
fig = autoplot(dBDI_log, size = .8, colour = 'black')
fig = fig +
theme_light() +
theme(plot.margin = ggplot2::margin(1, 1, 1, 1, "cm")) +
theme(text=element_text(size=30)) +
labs(x = "Month - Year") +
labs(y = "Transformed Baltic Dry Index (BADI)")
fig
SVAR modelΒΆ
A VAR model is preferred when analyzing multiple interdependent time series, as it allows for feedback and mutual interactions among variables, unlike the single-equation ADL model.
We tested stationarity since it ensures that relationships between variables remain stable over time. Non-stationary data can lead to spurious results, so differencing or a VECM formulation is needed if variables are cointegrated.